
import os
import glob
import pandas as pd

from matplotlib.gridspec import GridSpec

# Import cerebellum packages
import matplotlib.pyplot as plt
import SUITPy.flatmap as flatmap


bids_dir = '/data/projects/social_doors/'


os.chdir(bids_dir)

data_dir = os.path.join(bids_dir, 'derivatives','social_doors-nilearn')

# Define subject list
#subjs_scan_info = pd.read_csv(bids_dir+'/derivatives/mriqc/mriqc_summary_poor.csv')
#subjs_list = list(subjs_scan_info['subject'].unique())
#subjs_list.sort()
subjs_info = pd.read_csv(bids_dir+'derivatives/participants_good.tsv', sep='\t', index_col=0)
subjs_info = subjs_info.rename(columns={'participant_id': 'subject_label'})

# Remove participants with bad data
subjs_info = subjs_info[subjs_info['subject_label'].str.contains('sub-3880')==False]
subjs_info = subjs_info[subjs_info['subject_label'].str.contains('sub-4069')==False]


subjs_info_kids = subjs_info[subjs_info['group']=='kid']
subjs_info_colg = subjs_info[subjs_info['group']=='college']

subjs_list = subjs_info['subject_label'].to_list()
subjs_list_kids = subjs_info_kids['subject_label'].to_list()
subjs_list_colg = subjs_info_colg['subject_label'].to_list()

print('Found '+str(len(subjs_list_kids))+' adolescent subjects')
print('Found '+str(len(subjs_list_colg))+' college subjects')

Found 32 adolescent subjects
Found 29 college subjects


# Filter for releveant data
subjs_info_fltr = subjs_info[['subject_label', 'age']]
subjs_info_kids_fltr = subjs_info_kids[['subject_label','age']]
subjs_info_colg_fltr = subjs_info_colg[['subject_label','age']]


from nilearn.glm.second_level import make_second_level_design_matrix

design_matrix = make_second_level_design_matrix(subjs_list_kids, subjs_info_kids_fltr)

/opt/anaconda3/lib/python3.7/site-packages/nilearn/glm/__init__.py:56: FutureWarning: The nilearn.glm module is experimental. It may change in any future release of Nilearn.
  'It may change in any future release of Nilearn.', FutureWarning)
/opt/anaconda3/lib/python3.7/site-packages/nilearn/glm/first_level/design_matrix.py:475: UserWarning: Attention: Design matrix is singular. Aberrant estimates are expected.
  warn('Attention: Design matrix is singular. Aberrant estimates '


from nilearn.plotting import plot_design_matrix

plot_design_matrix(design_matrix)

<AxesSubplot:label='conditions', ylabel='scan number'>


from nilearn.glm.second_level import SecondLevelModel

mni_mask = bids_dir+"derivatives/fmriprep/sub-010/anat/sub-010_space-MNI152NLin2009cAsym_label-GM_probseg_bin.nii.gz"


contrast = 'valence_x_outcome'
task = 'social'
group = 'kids'

temp_file_list = []
for subj in subjs_list_kids:
    temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
    temp_file_list.append(temp_file[0])
    

model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=6.0)
model.fit(temp_file_list, design_matrix=design_matrix)

SecondLevelModel(mask_img='/data/projects/social_doors/derivatives/fmriprep/sub-010/anat/sub-010_space-MNI152NLin2009cAsym_label-GM_probseg_bin.nii.gz',
                 smoothing_fwhm=6.0)


z_map = model.compute_contrast('intercept', output_type='z_score')

z_map.to_filename(os.path.join(data_dir,'group',
                               'zmap_'+group+'_'+task+'_'+contrast+'.nii.gz'))


from nilearn.plotting import plot_stat_map

plot_stat_map(z_map, threshold=2, cut_coords=range(-65,66,10), display_mode='x')

/opt/anaconda3/lib/python3.7/site-packages/nilearn/plotting/img_plotting.py:300: FutureWarning: Default resolution of the MNI template will change from 2mm to 1mm in version 0.10.0
  anat_img = load_mni152_template()

<nilearn.plotting.displays._slicers.XSlicer at 0x7f1b4c0b8390>


from nilearn.glm import threshold_stats_img

_, threshold = threshold_stats_img(z_map, alpha=0.001, height_control='fpr')

plot_stat_map(z_map, threshold=threshold, cut_coords=range(-65,66,10), display_mode='x')

<nilearn.plotting.displays._slicers.XSlicer at 0x7f1b40b60390>


design_matrix = make_second_level_design_matrix(subjs_list_kids, subjs_info_kids_fltr)

#contrasts = ['valence_x_outcome',
#             'positive','positive_win','positive_loss',
#             'negative','negative_win','negative_loss']
contrasts=['all_winVlos', 'positive_winVlos']
contrasts_2nd = ['age','intercept']
tasks = ['mdoors','social']
group = 'kids'
alpha = 0.001


threshold_data = {}

for task in tasks:
    for contrast in contrasts:
        
        temp_file_list = []
        for subj in subjs_list_kids:
            temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
            temp_file_list.append(temp_file[0])
        temp_file_list.sort()
        
        print('Calculating '+task+' group '+contrast+' contrast')
        
        model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=8.0)
        model.fit(temp_file_list, design_matrix=design_matrix)
        
        for contrast_2nd in contrasts_2nd:
            z_map = model.compute_contrast(contrast_2nd, output_type='z_score')

            z_map.to_filename(os.path.join(data_dir,'group',
                                       'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_unc.nii.gz'))
            
            # Multiple Comparisons Correction
            z_map_thresh, threshold = threshold_stats_img(z_map, alpha=alpha, height_control='fpr')
            z_map_thresh.to_filename(os.path.join(data_dir,'group',
                                       'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr.nii.gz'))
            
            threshold_data['zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr_'+str(alpha)] = threshold

Calculating mdoors group all_winVlos contrast
Calculating mdoors group positive_winVlos contrast
Calculating social group all_winVlos contrast
Calculating social group positive_winVlos contrast


def plot_stat_map_cb(filename, threshold):
    # Find stat map
    stat_filename = os.path.join(data_dir,'group',filename+'.nii.gz')

    # Import stat map as a cerebellum flatmap
    funcdata = flatmap.vol_to_surf(stat_filename,
                                  space='SUIT')

    # Set figure specs
    fig = plt.figure(figsize=(15, 4))
    gs = GridSpec(2, 3)
    ax_img1 = plt.subplot(gs[0, :2])
    ax_img2 = plt.subplot(gs[1, :2])

    # Plot stat maps
    plot_stat_map(stat_filename, threshold=threshold, cut_coords=range(-65,0,10), display_mode='x', axes=ax_img1,
                 annotate=False, title='Left hemisphere')
    plot_stat_map(stat_filename, threshold=threshold, cut_coords=range(66,5,-10), display_mode='x', axes=ax_img2,
                 annotate=False, title='Right hemisphere')
    ax_joint = plt.subplot(gs[:, 2:])
    ax_joint.set(title='Cerebellum flatmap')
    flatmap.plot(data=funcdata, cmap='hot',
                threshold=[threshold, threshold+5],
                colorbar=True,
                render='matplotlib', new_figure=False)


filename = 'zmap_kids_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_kids_mdoors_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_kids_social_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_kids_social_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


design_matrix = make_second_level_design_matrix(subjs_list_colg, subjs_info_colg_fltr)


contrasts=['all_winVlos', 'positive_winVlos']
contrasts_2nd = ['age','intercept']
tasks = ['mdoors','social']
group = 'colg'
alpha = 0.001


threshold_data = {}

for task in tasks:
    for contrast in contrasts:
        
        temp_file_list = []
        for subj in subjs_list_colg:
            temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
            temp_file_list.append(temp_file[0])
        temp_file_list.sort()
        
        print('Calculating '+task+' group '+contrast+' contrast')
        
        model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=8.0)
        model.fit(temp_file_list, design_matrix=design_matrix)
        
        for contrast_2nd in contrasts_2nd:
            z_map = model.compute_contrast(contrast_2nd, output_type='z_score')

            z_map.to_filename(os.path.join(data_dir,'group',
                                       'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_unc.nii.gz'))
            
            # Multiple Comparisons Correction
            z_map_thresh, threshold = threshold_stats_img(z_map, alpha=alpha, height_control='fpr')
            z_map_thresh.to_filename(os.path.join(data_dir,'group',
                                       'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr.nii.gz'))
            
            threshold_data['zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr_'+str(alpha)] = threshold

/opt/anaconda3/lib/python3.7/site-packages/nilearn/glm/first_level/design_matrix.py:475: UserWarning: Attention: Design matrix is singular. Aberrant estimates are expected.
  warn('Attention: Design matrix is singular. Aberrant estimates '

Calculating mdoors group all_winVlos contrast
Calculating mdoors group positive_winVlos contrast
Calculating social group all_winVlos contrast
Calculating social group positive_winVlos contrast


filename = 'zmap_colg_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_colg_mdoors_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_colg_social_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_colg_social_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


subjs_info_num = subjs_info.copy()
subjs_info_num = subjs_info_num.replace({'sex': {'F': 0, 'M': 1},
                                         'group': {'college': -1, 'kid': 1}})
subjs_info_num = subjs_info_num.drop(columns=['age', 'sex'])

design_matrix = make_second_level_design_matrix(subjs_list, subjs_info_num)

plot_design_matrix(design_matrix)

<AxesSubplot:label='conditions', ylabel='scan number'>


#design_matrix = make_second_level_design_matrix(subjs_list, subjs_info_fltr)

#contrasts = ['valence_x_outcome',
#             'positive','positive_win','positive_loss',
#             'negative','negative_win','negative_loss']
contrasts=['all_winVlos', 'positive_winVlos']
contrasts_2nd = ['intercept','group']
tasks = ['mdoors','social']
group = 'all'
alpha = 0.001

threshold_all = {}

#temp_file_list = []
#for subj in subjs_list_kids:
#    temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii.gz'))
#    temp_file_list.append(temp_file[0])
for task in tasks:
    for contrast in contrasts:
        
        temp_file_list = []
        for subj in subjs_list:
            temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
            temp_file_list.append(temp_file[0])
        temp_file_list.sort()
        
        print('Calculating '+task+' group '+contrast+' contrast')
        
        model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=8.0)
        model.fit(temp_file_list, design_matrix=design_matrix)
        
        for contrast_2nd in contrasts_2nd:
            z_map = model.compute_contrast(contrast_2nd, output_type='z_score')

            z_map.to_filename(os.path.join(data_dir,'group',
                                       'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_unc.nii.gz'))
            
            # Multiple Comparisons Correction
            z_map_thresh, threshold = threshold_stats_img(z_map, alpha=alpha, height_control='fpr')
            z_map_thresh.to_filename(os.path.join(data_dir,'group',
                                       'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr.nii.gz'))
            
            threshold_data['zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr_'+str(alpha)] = threshold

Calculating mdoors group all_winVlos contrast
Calculating mdoors group positive_winVlos contrast
Calculating social group all_winVlos contrast
Calculating social group positive_winVlos contrast


filename = 'zmap_all_social_positive_winVlos_group_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_all_social_all_winVlos_group_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_all_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_all_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_all_social_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_all_mdoors_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


filename = 'zmap_all_social_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]

plot_stat_map_cb(filename+'unc', threshold)


behav_data = pd.read_excel(bids_dir+'derivatives/behavioral/aggregate_full_data_SPSS_020419.xlsx',
                          sheet_name='Sheet1')

behav_data.head()


behav_data_fltr = behav_data[['Subject','ch_totanx', 'ch_cdi_total']]

behav_data_fltr['Subject'] = behav_data_fltr['Subject'].str.replace('s','sub-')

behav_data_fltr = behav_data_fltr.rename(columns={'Subject':'subject_label'})

behav_data_fltr.head()

/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  This is separate from the ipykernel package so we can avoid doing imports until


# Add behavioral data to design matrix

subjs_info_kids_behav = subjs_info_kids_fltr.merge(behav_data_fltr, on='subject_label')


design_matrix = make_second_level_design_matrix(subjs_list_kids, subjs_info_kids_behav)

plot_design_matrix(design_matrix)

<AxesSubplot:label='conditions', ylabel='scan number'>

	Order	Subject	ch_totanx	ch_socialanx	p_totanx	p_socialanx	bfne	ch_cdi_total	ch_cdi_interpersonal	ch_cdi_anhedonia	...	FB_VS_6mm_R_Like_Lose_NZMean_1	FB_VS_6mm_R_Like_Lose_NZCount_1	FB_VS_6mm_R_Like_Lose_NZMin_1	FB_VS_6mm_R_Like_Lose_NZMax_1	FB_VS_6mm_R_Like_Win_Mean_1	FB_VS_6mm_R_Like_Win_NZMean_1	FB_VS_6mm_R_Like_Win_NZCount_1	FB_VS_6mm_R_Like_Win_NZMin_1	FB_VS_6mm_R_Like_Win_NZMax_1	FB_VS6mm_L_Correct_DislikeVsLike
0	1	s010	51	14	7	1	33	11.0	0.0	5.0	...	0.029036	123.0	-0.296210	0.647286	0.256033	0.256033	123.0	-0.022722	0.752379	-0.127002
1	1	s011	12	3	4	1	31	8.0	0.0	3.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	1	s013	8	4	1	1	13	5.0	0.0	2.0	...	-0.308424	123.0	-0.554069	0.030482	0.093248	0.093248	123.0	-0.183241	1.003483	-0.167305
3	1	s015	18	5	1	0	21	NaN	NaN	NaN	...	-0.798996	123.0	-1.841710	-0.031699	0.004660	0.004660	123.0	-1.020790	0.809193	0.485397
4	2	s024	21	10	0	0	28	1.0	0.0	0.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

	subject_label	ch_totanx	ch_cdi_total
0	sub-010	51	11.0
1	sub-011	12	8.0
2	sub-013	8	5.0
3	sub-015	18	NaN
4	sub-024	21	1.0

Table of Contents

Adolescents¶

Second Level Design Matrix¶

Valence x Outcome Sanity Check¶

Adolescent Cohort¶

Monetary Task¶

Positive win > Positive loss¶

All wins > All losses¶

Positive win > Positive loss¶

All wins > All losses¶

Young Adult Cohort¶

Monetary Task¶

Positive win > Positive loss¶

All wins > All losses¶

Positive win > Positive loss¶

All wins > All losses¶

All Subjects¶

Group differences¶

Monetary Task¶

Positive win > Positive loss¶

All wins > All losses¶

Accounting for Age¶

Behavioral Analysis¶

Table of Contents

Adolescents¶

Second Level Design Matrix¶

Valence x Outcome Sanity Check¶

Adolescent Cohort¶

Monetary Task¶

Positive win > Positive loss¶

All wins > All losses¶

Social Task¶

Positive win > Positive loss¶

All wins > All losses¶

Young Adult Cohort¶

Monetary Task¶

Positive win > Positive loss¶

All wins > All losses¶

Social Task¶

Positive win > Positive loss¶

All wins > All losses¶

All Subjects¶

Group differences¶

Monetary Task¶

Social Task¶

Positive win > Positive loss¶

All wins > All losses¶

Accounting for Age¶

Behavioral Analysis¶